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Abstract 

We  study  the  problem  of  estimating  the  log  spectrum  of  a  stationary  Gaussian 
time  series  by  thresholding  the  empirical  wavelet  coefficients.  We  propose  the  use  of 
thresholds  depending  on  sample  size  n,  wavelet  basis  ^  and  resolution  level  j.  At 
fine  resolution  levels  (j=  1,2, ...),  we  propose 

loS 

where  {qj}  are  level-dependent  constants  and  at  coarse  levels  (j  >•  1) 

The  purpose  of  this  thresholding  level  is  to  make  the  reconstructed  log-spectrum  as 
nearly  noise-free  as  possible.  In  addition  to  being  pleasant  from  a  visual  point  of  view, 
the  noise-free  character  leads  to  attractive  theoretical  properties  over  a  wide  range  of 
smoothness  aissumptions.  Previous  proposals  set  much  smaller  thresholds  and  did  not 
enjoy  these  properties. 
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1  Introduction 


Suppose  we  want  to  study  a  time  series  from  an  observed  segment  of  the  series  X\,...,Xn- 
There  are  two  aspects  to  the  study  of  time  series  —  analysis  and  modelling.  The  aim  of 
analysis  is  to  understand  the  nature  of  the  process.  The  main  reason  for  modelling  a  time 
series  is  to  enable  forecasts  of  future  values  of  the  process.  No  forecasting  caa  be  made 
before  we  understand  the  salient  features  of  the  process  being  characterized. 

The  analysis  of  a  time  series  can  be  done  either  in  the  time  domain  or  in  the  frequency 
domain.  In  the  time  domain,  attention  is  focused  on  the  relationship  between  observations 
at  different  points  in  time,  e.g.  ARMA  models;  while  in  the  frequency  domain  it  is  cyclical 
movements  which  are  studied  and  this  can  be  done  by  studying  the  spectral  density  h. 
The  two  forms  of  analysis  are  complementary  rather  than  competitive.  They  give  different 
insights  into  the  nature  of  the  process. 

Our  goal  in  this  paper  is  to  estimate  h  from  the  data  Xi,  ...,Xjv.  Here 

1  °° 

^  7(5)cos(sa;). 

5=2—00 


The  techniques  currently  used  in  spectral  analysis  are: 

•  Window  methods  (also  known  as  kernel  methods!.  §6.2.3  of  Priestley  (1981).  A  spec¬ 
tral  window  (also  called  weight  function)  is  used  to  smooth  the  periodogram,  or 
equivalently  a  lag  window  (also  known  as  weight  sequences,  Priestley  (1981),  pp. 
437,  covariance  window,  Bentkus  and  Susinskas  (1982))  W(-)  is  applied  to  the  sample 
auto-covariances.  The  resulting  estimate  is: 

=  ^  S  W^(s)7(5)cos(sa;) 

|s|<iV 

Typical  choices  of  windows  include  the  Fejer  window:  W(s)  =  -f[|3|<Af]  for  some 
1  <  M  <  N ,  for  example  M  =  y/N ;  Bartlett  window:  W{s)  =  (1  -  |s|/M)+;  Daniell 
window,  W(s)  =  sin(7rs/M)/(Ts/M);  see  §6.2.3  of  Priestley  (1981)  for  more  choices 
of  windows.  Bentkus  and  Susinskas  (1982)  has  studied  optimal  L2  convergence  rate 
for  the  window  method  over  certain  classes  of  smooth  spectra. 


•  Mixtures  of  the  above,  e.g.  the  prewhitening  technique,  §7.4.1  of  Priestley  (1981). 

•  Maximum  entropy  (ME)  methods,  introduced  by  Burg  (1967,  1972).  This  method  is 
essentially  equivalent  to  AR  method. 

•  ARMA  approximation:  similar  idea  as  AR  approximation. 

They  perform  well  in  many  cases  and  have  been  used  frequently  in  practice.  They  have 
both  advantages  and  disadvantages.  For  example,  window  methods  have  computational 
advantages  but  they  may  perform  poorly  in  cases  of  high  dynamic  range;  AR  approxima¬ 
tion  is  suitable  only  for  relatively  smooth  spectra,  Tukey  (1978);  ARMA  procedures  are 
almost  always  effective  and  furnishes  attractive  interpretations  for  the  derived  model;  but 
it  is  computationally  unrealistic  to  fit  high-order  ARMA  models,  and  rigorous  theory  of 
asymptotic  properties  is  lacking. 

In  this  paper,  we  will  develop  a  technique,  based  on  wavelet  decomposition  of  the 
periodogram  and  reconstruction  of  the  spectrum;  our  new  approach  avoids  shortcomings 
of  earlier  methods.  Our  new  approach  is  computationally  efficient;  it  can  estimate  spectra 
which  are  nonsmooth  at  a  near-optimal  rate;  and  it  can  be  rigorously  analyzed. 

It  is  well  known  that  the  asymptotic  variances  of  the  empirical  wavelet  coefficients 
are  functions  of  the  unknown  spectrum,  and  that  non-linear  Wavelet  Shrinkage  procedure 
below  heavily  depends  on  the  asymptotic  variances.  Since  the  asymptotic  variances  are 
“proportional”  to  the  coefficients,  i.e.,  the  bigger  the  coefficients  at  location  {j,k),  the 
higher  the  noise  level  at  this  location.  This  makes  estimation  more  difficult  in  practice.  In 
order  to  overcome  the  heteroscedasticity,  we  therefore  consider  estimating  the  log-spectrum. 

Suppose  we  observe  a  segment  Xo,Xi,  ...,X2n-i  of  a  Gaussian  time  series  with  mean 
zero  and  spectral  density  h,  where  for  simplicity  we  assume  that  n  =  2"*  is  dyadic.  In 
this  paper  we  discuss  the  following  4-step  wavelet  shrinkage  procedure  for  estimation  of  log 
spectrum,  g  =  logh: 


[1]  Calculate  the  log-periodogram 


Zi  =  log/(u;;),  /  =  0,...,n-l. 


where  ui  =  l/2n  and 


liu) 


1 

4ra7r 


2n-l 


t=o 


(1) 


[2]  Take  a  standard  periodic  wavelet  transform  of  (z;)  to  get  the  empirical  wavelet  coef¬ 
ficients  {yj,k}j=l,2,...m—l,k=0,...,n/V  —  l' 

[3]  Apply  the  soft  threshold 

St{x)  =  sgn(x)(|a:l  -  t)+  (2) 
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to  the  empirical  wavelet  coefficients  {yj,k}^  with  level-dependent  thresholds 
according  to  formula 

tj,n  =  -^V^ogny  ajlogn 


(3) 


where  (oi , aio)  =  ( 1.29, 1.09, 0.92, 0.77, 0.65, 0.54, 0.46, 0.39, 0.32, 0.27)  for  the  com¬ 
monly  used  compactly  supported  orthogonal  wavelet  bases  {coiflets,  doublets  and 
symmlets,  Chapters  6  and  8  of  Daubechies,  1992)  and  some  selected  n. 


[4]  Invert  the  wavelet  transform,  producing  an  estimate  (5*)  of  the  log-spectrum  at  the 
Fourier  frequencies  ui. 


This  procedure  falls  in  the  general  category  of  wavelet  shrinkage  estimates  for  noisy 
data.  Donoho  and  Johnstone  (1992a,b,c),  have  discussed  the  application  of  such  thresholded 
wavelet  transforms  for  recovering  curves  from  noisy  data,  where  the  noise  is  assumed  to 
have  a  Gaussian  distribution,  and  they  have  proposed  the  level-independent  threshold  tn  = 
cri/21og(n),  where  a  is  the  noise  level.  Such  methods  have  also  been  developed  in  density 
estimation:  Johnstone,  Kerkyacharian  and  Picard  (1992)  propose  the  use  of  level- dependent 
thresholds  ^/J.  In  general  such  wavelet  shrinkage  methods  have  a  number  of  theoretical 
advantages,  including  near-optimal  mean-squared  error  ajid  near-ideal  spatial  adaptation. 

Recently,  the  author  and  others  have  studied  the  possibility  of  using  wavelet  shrinkage 
to  de-noise  periodogram  data.  Application  of  wavelet  shrinkage  to  the  log-periodogram  is 
particularly  attractive,  Moulin  (1993a,  b),  since  the  logarithm  is  the  variance  stabilizing 
transformation  for  the  periodogram,  Wahba  (1980). 

It  is  easy  to  apply  wavelet  shrinkage  software  designed  for  Gaussian  noise  to  the  log  pe¬ 
riodogram,  and  in  some  cases  acceptable  reconstructions  have  been  obtained.  However,  this 
will  not  always  be  the  case.  The  noise  in  the  wavelet  coefficients  of  the  log-periodogram  has 
a  non- Gaussian  character;  it  reaches  high  levels  somewhat  more  frequently  than  Gaussian 
theory  would  predict.  Consequently,  thresholds  set  based  on  Gaussian  theory  will  not  be 
high  enough  to  completely  suppress  noise  in  the  coefficients.  The  thresholding  we  propose 
is  based  on  a  careful  analysis  of  the  non-Gaussian  character,  and  is  somewhat  larger  than 
Gaussian  theory  would  predict. 

As  an  example  of  the  difference  between  our  method  and  other  proposals,  we  present 
in  Figure  1  a  display  for  an  Ai2(24)  time  series:  (a)  its  log  spectrum;  (b)  log  periodogram; 
and  wavelet  reconstructions  based  on  (c)  our  proposal  and  (d)  based  on  Gaussian  theory.  It 
is  visually  evident  that  the  spectral  estimate  based  on  Gaussian  theory  has  spurious  noise 
spikes,  which  are  unrelated  to  the  true  underlying  spectrum.  More  plots  for  this  example, 
along  with  other  examples,  for  different  sample  sizes  and  different  wavelet  bases,  are  given 
at  the  end  of  the  paper. 

The  pattern  visually  evident  in  this  example  is  confirmed  by  two  theorems  which  we 
prove  in  this  paper. 
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Theorem  1  Under  the  Wahba  approximation  for  log  periodograms  (see  section  2.2  below), 
05  n  — f  00, 


P{U[sup  \yi,k  -  Eyj,k\  >  tj,n]}  0.  (4) 

In  words,  the  thresholds  are  set  high  enough  so  that  the  noise  does  not  surpass  them. 

Theorem  2  Under  the  Wahba  approximation  for  log  periodograms  and  a  compactly  sup¬ 
ported  wavelet  basis,  at  the  finest  level  j  =  l,  as  n  oo, 

P{sup  \yi,k  -  Eyi,k\  >  \/21ogn}  -h-  1.  (5) 

k 

Similar  results  hold  at  levels  j  =  2, 3, .... 

In  words:  noise  spikes  are  almost  surely  to  be  absent  for  our  proposal  and  almost  surely 
to  be  present  in  a  proposal  based  on  Gaussian  noise. 

The  significance  of  these  noise  spikes  and  the  need  to  de-noise  is  more  than  cosmetic.  In 
a  recent  article,  Donoho,  Johnstone,  Kerkyacharian  and  Picard  (1992)  study  the  Gaussian 
white  noise  model,  and  show  that  by  exploiting  a  noise-free  property  like  (4)  one  can 
show  that  the  estimator  attains  near-optimal  reconstruction  in  a  wide  variety  of  norms 
simultaneously  over  a  wide  variety  of  smoothness  assumptions.  Their  arguments  are  general 
and  abstract  and  transfer  to  the  present  setting  after  one  has  established  the  de-noising 
property  (4) 

If  noise-free  property  like  (4)  is  not  maintained,  then  an  inspection  of  their  arguments 
will  show  that  the  resulting  estimate  does  not  achieve  near  optimal  reconstruction  in  norms 
which  measure  smoothness  of  the  reconstruction,  although  it  may  well  be  true  that  the  error 
in  £2  norm,  which  does  not  measure  smoothness,  is  stiU  acceptable. 

2  The  Choice  of  Thresholds 

Our  choice  of  the  thresholds  (tj.n)  depends  on  the  tail  behavior  of  the  empirical  coefficients 

{yj,k}-  Let 

Pj{t)  =  m^P{\yj,k  -  Eyj,k\  >  t}  (6) 

be  the  tail  probability  of  where  wavelet  coefficient  yj^k  is  a  standardized  linear  combi¬ 
nation  of  the  log  periodogram  ordinates.  Then  we  aim  to  establish  that 

n  Pj(,lj,n}  ^  0)  n.  ^  00,  ( I ) 

uniformly  for  j  =  1, 2, ...,  m  -  1  as  this  implies  (4)  in  Theorem  1  (see  (15).  In  this  section, 
we  would  like  to  give  two  heuristic  arguments,  upon  which  our  proposed  thresholds  (3)  are 
based. 
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2.1  Normal  Approximation 

For  normally  distributed  wavelet  coefficients 

Wj^k  —  ^i,k  +  £j,k 

where  Sj^k  ~  A’(0,o’^)  are  iid,  the  threshold 

=  tn=  <T\/21ogn 

enjoys  a  variety  of  nice  theoretical  properties,  Donoho  and  Johnstone  (1992a,  b,  c).  Under 
certain  regularity  conditions,  for  j  close  to  m  =  [log2  n]  (technically  j  -+  oo  as  m  — *•  oo), 
as  n  — s-  oo,  individual  empirical  wavelet  coefficients  will  be  asymptotically  normal: 

Wj^k  -  Ewj^k  ~  A(0,  TT^/e) 

this  follows  by  applying  results  of  Taniguchi  (1979,  1980).  Moulin  (1993a,  b)  makes  similar 

argument.  So  when  the  sample  size  n  is  large,  for  coarse  levels  (j  close  to  m), 

are  approximately  normally  distributed  with  the  same  variance  tt^/S.  One  may  argue  that 

«(’"\/l°S(”)/^)  S  "•  >  J  ’  (8) 

Therefore  (7)  holds  and  this  leads  to  the  first  part  of  our  thresholds. 

Here  “regularity”  refers  to  the  length  of  memory  and  the  asymptotic  normality  for  the 
empirical  coefficients  not  only  depends  on  n  (large)  and  j  (large),  but  also  the  regularity  of 
the  time  series.  In  general,  the  longer  the  memory,  the  weaker  the  asymptotic  normality. 

2.2  Wahba  Approximation 

For  small  j,  the  normal  approximation  deteriorates.  Under  a  compactly  supported  wavelet 
basis  (chapters  6  and  8  of  Daubechies,  1992),  for  fixed  j,  the  coefficients  of  the  finest  level 
yj^s  are  linear  combinations  of  fixed  number  of  {log/(u;;)}.  For  example,  with  the  Haar 
wavelet,  yi^k  is  just  the  difference  of  two  adjacent  log  J(a;;)’s. 

Consider  the  following  non-Gaussian  additive  noise  model: 

zi  =  gi  +  ei  (9) 

I  —  1,2,  ...,n  -  1,  where  (gi)  is  the  object  to  be  estimated  (e.g.  log  spectrum  at  frequency 
a;;)*and 


£,•  =  log(77,y2)  +  7  (10) 

where  {r/,-,  i  —  1, 2, ...,  n  -  1}  are  independently  xi  =  exp(l/2)  distributed  and  7  =  0.57721 
is  the  Euler-Mascheroni  constant.  It  can  be  shown  that  Esi  =  0  and  var(£,)  =  -^/6. 
Wahba  (1980)  proposed  this  as  a  model  for  the  log  periodogram,  i.e.,  zi  w  logI{u!i)  and 
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gi  RJ  log/(u;;)  -  7.  Note  that  I{uq)  ~  Xi-  Since  the  influence  of  this  term  is  negligible  for 
n  lajge,  we  will  ignore  this  term  in  our  discussion. 

For  circulant  time  series,  see  §4.3  of  Harvey  (1989),  log/(a;/)  follows  model  (9)-(10) 
exactly.  For  general  stationary  processes,  the  above  model  is  only  asymptotically  true  (see 
Theorem  5.2.6  of  Brillinger,  1981).  The  exact  distribution  of  log  periodogram  can  be  found 
in  Wittwer  (1986). 


Now  we  will  show  that  under  the  model  (9)-(10),  our  advertised  thresholds  are  indeed 
adequate  for  the  finer  levels.  Note  that  under  model  (9)-(10),  each  yj^k  is  a  standardized 
lineax  combination  of  Zi’s. 

Let  us  use  the  very  finest  level,  i  =  1,  as  an  example.  Similar  results  hold  for  levels 
j  =  2, 3, ....  Suppose  a  compactly  supported  wavelet  basis  with  L  coefficients  in  the  dilation 
equation.  For  example,  for  the  Haar  wavelet,  X  =  2;  D4  wavelet,  X  =  4,  etc.,  more  examples 
can  be  found  in  Chapters  6  and  8  of  Daubechies  (1992).  Then 

L 

yi,k  =  ^  aiZi+2k-3 
/=! 

where  s  is  some  shift  parameter  for  computational  purpose  (e.g.  s  =  i/2)  and  — 
0,  12®?  =  Inequalities  below  leads  to  the  threshold  (3). 

Let  a  =  maxi</<i  |a/|,  from  the  proofs  of  the  Theorems  (see  (16)  and  (22)  of  Appendix), 
we  can  show  that  for  t  large, 

0.25€“^‘^“  <  p\{t)  <  {etj .  (11) 


Hence,  for  n  large,  there  is  a  constant  A  so  that 


<Pi(alogn)< 


„2a/a  ^  „a/a  '  (^2) 

Hence  thresholds  at  the  finest  level  must  have  the  form:  alogn,  for  some  constant  a  >  a, 
and  then 


npi(Q:log(2Ti))  0,  n  oo.  (13) 

Note  that  the  wavelet  filters  obey  ®?  =  !>  which  implies  a  <  1.  Hence,  for  the  finest 
level,  taking  a  =  1  always  guarantees  (13)  and  (7)  follows  for  our  proposal 
In  general,  there  are  Lj  non-zero  constants  a,j  so  that 

yj,k  —  (14) 

where  ©  in  subscript  is  interpreted  modulo  n.  Let  aj  =  majq  |a,j|,  then  we  have 

lj  <  2\L  -  1)  -  i  -I-  2  and  a,  <  a2-(■'-^)/^ 

where  L  is  the  length  of  the  filter  and  a  =  supj>i  <  1.534.  In  Tables  1  and  2  of 

Appendix,  we  list  some  Lj's  and  2f-^"^)/^aj’s  for  the  commonly  used  wavelets  (apparently 
they  converge  to  sup^,  |^(x)|).  This  leads  to  the  threshold  (3). 
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3  Numerical  Examples 

In  this  section,  we  study  four  time  series  and  compare  our  proposed  method  with  the  one 
based  on  Gaussian  theory. 

Figure  1  contains  an  AB,(24)  signal,  with  coefficients: 

-2.5216281  4.7715359  -7.9199915  11.9769211  -16.0778828 

20.6343346  -25.0531521  28.8738136  -31.8046265  34.0071373 

-34.7700272  34.3151321  -32.7861099  30.2861233  -26.7109356 

22.8838310  -18.7432098  14.5717688  -10.7177744  7.5322194 

-4.7226319  2.6807923  -1.3391306  0.5167125 

Figure  2  contains  a  white  noise  sigani,  the  true  log  spectral  density  in  this  case  is  giu)  = 
log(cr^/27r),  a  constant. 

Figure  3  contains  an  MA(15001)  time  series  with  coefficients 

sin(xn/2)  ,  ^ 

fill  =  7r/4,  a„+i  = - ,  n  =  1, 2, ...,  15000. 

n 

AR  method  is  compared  in  this  example. 

Figure  4  contains  the  Sunspots  signal,  WaveShrinh  estimate  with  comparisons  with  an  AR 
estimate  from  SPLUS  function  spec.ar.  An  AR(2)  model,  suggested  by  Priestley  (1981), 
pp.  882,  is  also  plotted. 
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(C)  (0) 

Figure  1:  (a),  log  spectrum  of  an  AR(24)  time  series,  (b).  log  periodogram,  (c).  WaveShrink 
estimate  based  the  new  thresholding  scheme,  and  (d).  WaveShrink  estimate  based  on  Gaussian 
theory. 
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(0  W) 


Figure  2:  (a),  log  spectrum  of  a  white  noise  signal,  (b).  log  periodogram,  (c).  WaveShrink  estimate 
based  the  new  thresholding  scheme,  and  (d).  WaveShrink  estimate  based  on  Gaussian  theory. 
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Figure  3:  (a),  log  spectrum  of  a  long  MA  series,  (b).  log  periodogram,  (c).  WaveShrink  estimate 
based  the  new  thresholding  scheme,  and  (d).  AR  estimate  from  SPLUS  function  spec.ar. 
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Figure  4:  (a),  sunspots  data,  (b).  log  periodogram,  (c).  WaveShrink  estimate  (solid  line)  and 
AR(2)  estimate  (dotted  line),  and  (d).  AR  (of  order  27)  estimate. 
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4  Appendix 

4.1  Tables 


IB 

LI 

HI 

L2 

H2 

L3 

H3 

L4 

H4 

L5 

H5 

L6 

H6 

C6 

0.85 

0.85 

0.95 

1.10 

1.02 

1.28 

1.07 

1.40 

1.10 

1.48 

C12 

0.81 

0.81 

0.85 

1.01 

0.86 

1.08 

0.87 

1.10 

0.87 

1.11 

0.87 

1.12 

CIS 

0.79 

0.79 

0.81 

0.96 

0.81 

0.99 

0.81 

1.00 

0.81 

1.00 

0.81 

1.00 

C24 

0.78 

0.78 

0.79 

0.93 

0.79 

0.94 

0.79 

0.95 

0.79 

0.95 

0.79 

0.95 

C30 

0.77 

0.77 

0.78 

0.91 

0.78 

0.91 

0.78 

0.92 

0.78 

0.92 

0.78 

0.92 

D2 

0.71 

0.71 

0.71 

0.71 

0.71 

0.71 

0.71 

0.71 

D4 

0.84 

0.84 

0.95 

1.08 

0.99 

1.23 

1.00 

1.28 

1.00 

1.29 

0.99 

1.28 

D6 

0.81 

0.81 

0.83 

0.98 

0.86 

1.03 

0.89 

1.11 

0.90 

1.15 

0.90 

1.18 

D8 

0.71 

0.71 

0.79 

0.77 

0.78 

0.95 

0.79 

0.96 

0.79 

0.96 

DIO 

0.72 

0.72 

0.81 

0.94 

0.74 

0.89 

0.75 

0.83 

0.76 

0.84 

0.76 

0.84 

D12 

0.75 

0.75 

0.75 

0.92 

0.74 

0.80 

0.73 

0.80 

0.73 

0.8 

D14 

0.73 

0.73 

0.74 

0.79 

0.73 

0.82 

0.73 

0.79 

0.73 

0.79 

0.73 

0.79 

D16 

0.68 

0.68 

0.74 

0.82 

0.71 

0.82 

0.71 

0.79 

0.71 

0.78 

0.71 

0.78 

D18 

0.66 

0.66 

0.72 

0.87 

0.70 

0.74 

0.70 

0.75 

0.70 

0.76 

0.70 

0.76 

D20 

0.69 

0.69 

0.69 

0.75 

0.68 

0.73 

0.69 

0.72 

0.69 

0.73 

S8 

0.99 

0.85 

1.05 

0.85 

1.05 

0.86 

1.07 

SIO 

0.72 

0.72 

0.77 

0.77 

0.78 

0.89 

0.79 

0.93 

0.79 

0.94 

0.79 

0.95 

S12 

0.79 

0.79 

0.80 

0.95 

0.80 

0.97 

0.81 

0.97 

0.81 

1.00 

0.81 

1.00 

S14 

0.77 

0.77 

0.79 

0.92 

0.77 

0.92 

0.77 

0.89 

0.77 

0.89 

S16 

0.78 

0.78 

0.78 

0.91 

0.78 

0.92 

0.79 

0.93 

0.79 

0.95 

S18 

0.72 

0.72 

0.76 

0.76 

0.77 

0.89 

0.77 

0.77 

0.89 

0.77 

0.89 

S20 

0.77 

0.77 

0.77 

0.89 

0.77 

0.89 

0.78 

0.91 

0.78 

0.92 

0.78 

0.92 

Table  1:  Maximum  Coefficients,  Lj  low-pass  cascade  filter  at  level  j,  Hj  high-pass  cascade 
filter  at  level  j.  C  for  Coiflets;  D  for  Daublets,  Daubechies’  original  wavelets;  S  for  Symmlet. 
Daubechies’  near-symmetric  wavelets. 
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L=2 

L=4 

L=6 

00 

II 

L=10 

II 

L=14 

L=16  L=18 

L=20 

j=l 

2 

4 

6 

8 

10 

12 

14 

16 

18 

20 

j=2 

4 

10 

16 

22 

28 

34 

40 

46 

52 

58 

j=3 

8 

22 

36 

50 

64 

78 

92 

106 

120 

134 

j=4 

16 

46 

76 

106 

136 

166 

196 

226 

256 

286 

j=5 

32 

94 

156 

218 

280 

342 

404 

466 

528 

590 

j=6 

64 

190 

316 

442 

568 

694 

820 

946 

1072 

1198 

Table  2:  Number  of  Non-zero  Coefficients,  L  is  the  length  of  the  filter  and  j  =  level. 

4.2  Proofs  of  the  Theorems 
Proof  of  Theorem  1:  First  of  aU, 

m— 1 

P{U[sup  \yj,k-Eyj,k\ >  tj,n]}  <  X) 

j  ^  i=i  * 

m— 1 

l<j<m 

Under  model  (9)  and  (10),  from  (14),  coefficients  in  level  j  can  be  written  ais 

yj,k  ~  Eyj,k  —  XI 

t=i 

with  —  0  3,nd  ^l,j  —  1-  =  maXt 

The  moment  generating  function  for  e  ~  log(x2/2)  is 

M(u)  =  =  r(u  +  1) 

where  T  is  Euler’s  Gamma  function.  Let 

A“(t)  =  inf{-st  +  XI(logr(l  +  Uijs)  +  logr(l  -  a.js))} 

“  j=i 

It  can  be  shown  that  for  any  random  variable  X  with  EX  =  0, 

♦  P{X  >a)< 


for  all  a  >  0  and  t.  Then 

Pi(0  < 


and  therefore, 


max  p,(t)  <  sup  e^J^‘^  =  e.xp(  sup  ANt)) 

l<j<m  l<j<Tn 


(16) 
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Let  In  =  (log  n)^/®  and  consider 


^i  = 


a  log"  T  ■  c'  f 

,1/4  ^ 

_ 

/SVlogra  Xj  >  4 


For  the  thresholds 


we  have  Therefore, 


Li 


and 


sup  Pi(ti,n)  <  sup  pjit])  <  exp(  sup 

l<j<m  l<j<m 


sup  A*(t*)  <  sup  A-(t*)V  sup  A;(tp 

l<j<m  l<Lj<in  t„<L,<n 


=  .  V.  AK^v^°s^) 


l<I,j<r„  '  X-'’  ’  U<L,<n 


When  Xj  <  in  =  (log  and  ra  large,  we  have 


a  log  n  /7“  ^  - 

ri/4  > 


Consider  the  following  non-positive  function  H,  defined  on  R+  by 

E{t)  =  -t  4- 1  -}-log< 


(17) 


(18) 


Lemma  1  For  t  >  ajLj, 

a;(0  <  ii-ffCrr-) 

Uj  J 

and 

X  i  in  X  for  large  n  and  1  <  L  <in. 

The  proof  wiU  be  given  later.  Combining  (16)  and  Lemma  1  yields  the  second  part  of  (11). 
From  the  Lemma, 


* 


,.,alogn^  ^ 

sup,  ^ji-znr)  ^ 


ir 


^  ^,alogn, 
sup 

aL/ 


<  5'(-logn)  = --logn  (1  +  0(1)) 
a  a 


When  Lj  >  ^n,  we  have 

/?>/logn  < 

and 


(19) 


15 


Lemma  2  For  Sj  €  (0, 1)  and  0  <  t  <  ir'^Sj/Saj, 


A -(0  <  - 


x2/3  +  6^; ' 


In  particular,  for  n  large  and  Sj  =  x  <  (logn)“^/^^, 

A-(()<-^(l  +  »(l))  t<i3Lf. 

The  proof  will  be  given  later. 

From  the  Lemma, 


sup  hIj{l5\/[og  n)  <  — 

tn<L,<n 


3/3^  log  n 


(1  +  0(1)). 


Combining  (18),  (19)  and  (20),  we  have 

sup  A;(tj-)  < -(^/\^)logn(l  +  o(l)). 

l<J<m  «  ~ 

Therefore,  for  any  a>  a  and  /3  >  x/\/Z, 


n  sup  pj(tj,n)  <  nexp(  sup  A^(ij,„))  =  o(l). 

l<j<m  l<j<m 

and  this  completes  the  proof  of  Theorem  1. 


(20) 


Proof  of  Lemma  1:  First  of  all,  it  can  be  shown  that  for  0  <  x  <  1,  r(l  +  x)  <  r(l  —  z) 
and  r(x)  <  1/x.  (In  fact,  r(x)  =  x~^  flfei  +  1/^)^  <  1  +  when 

0  <  X  <  1.)  Therefore, 


A*(t)  <  inf  {-st  +  £logr(l  +  aijs)} 

aj 


i=l 


,  .  r  St 

<  Lj  mf  { - — 

■'  o<s<i  UjLj 


-log(l-s)} 


and  the  first  part  of  the  Lemma  can  be  easily  solved  from  the  last  expression. 

For  n  large,  such  that  Cn  =  aa“^logn  >  11  and  In  —  (logn)^/®  <  (Cn/ll)"*^^-  Let 
h{x)  =  xH{CnX~^^^),  then  for  1  <  x  <  f„,  h'{x)  is  increasing  in  x  and  h'iin)  <  0. 
Therefore,  h'{x)  <  0  for  1  <  x  <  and  this  means  h{x)  is  decreasing  in  x.  The  proof  is 
completed. 


Proof  of  Lemma  2:  Let  K{x)  =  logr(l  +  x),  then  for  x  >  -1,  by  Taylor  expansion, 

A"(x)  =  A(0)  +  /i:'(0)x  +  x^K"{ex)/2  =  -7X  +  x^K"i9x)/2,  (21) 
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where  7  =  0.5772....  is  the  Etiler  constant  and  9  =  9(x)  €  (0, 1).  Note  that  K'(x  -  1) 
is  the  famous  Psi  function.  By  the  properties  of  Psi  function  (Davis  (1935),  pp.  11), 
ii^"(a:)  =  E^=i(k  +  ®)~^  and  K"(0)  =  ttVO.  Hence,  for  any  S  €  (0, 1),  we  have  K"i-6)  < 
K"{0)  +  ZS  =  Ks  and  for  x  >  —S, 

K{x)  <  -71  +  x^Ks/2. 

Recall  that  aij  =  0  and  ajj  =  1, 

A*(t)  =  M{-st  +  '^Kiaijs)} 

i 

<  inf  {-St  +  +  aljs'^Ks/2)} 

i 

=  inf  {—st  +  s^Ks/2} 

\A<s,hi 

and  Lemma  2  follows  immediately. 

Proof  of  Theorem  2:  To  prove  Theorem  2,  we  need  a  lower  bound  for  Pj(t).  It  is  easy 
to  show  that  if  X,  Y  are  independent,  then  for  any  s,t  >  0, 

P{|X  +  y|  >  0  >  P{|Xj  <  s}P{|y|  >  t  +  s}. 

Suppose  |ai|  =  max  |a/|  =  a.  There  exists  to  >  0  such  that  for  t  >  to, 

^  1 

P{|2^a(log(j7;/2)|  <  t}  >  -. 

/=2  ^ 

Also  notice  that  for  0  <  i  <  1,  1  —  e~^  >  x/2.  Then 

L 

Pi(t)  =  P{|^a,log(77//2)|  >  t} 

;=i 

L 

>  P{lailog(7i/2)|  >  2t}P{|5]ojlog(7//2)|  <  t} 

1=2 

>  ^■P{Uog( 771/2)1  >  ^} 

>  ^^{log(7?i/2)  < 

=  i{l-exp(-e-«/»)}>ie-2‘/“.  (22) 

and  this  completes  the  proof  of  Lemma  1. 

From  (12)  we  can  see  that  for  any  constant  C  <  a  12, 

npi(Clogn) -H- 00,  n  — ►  00. 
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In  particular, 


npi{\/2logn)  =  npi(. 


logn 


logn)  -V  00, 


n 


00. 


There  exists  K  >  0  such  that  {pi.jk}  are  independent  (linear  combinations  over  disjoint 
intervals).  Therefore, 


P{sup  \yi,k  -  Eyi,k\  >  >72 logn}  > 
k 


P{sup  \yi,iK  -  Eyi^iKl  >  ^^logn} 

i 

n/2K 

1-  n  (1  -pi(\/2iogn))  1, 

t=0 


and  this  implies  the  Theorem  2. 

References 

[1]  Bentkus,  R.  Yu.  and  Susinskas,  Ju.  V.  (1982)  On  Optimal  Statistical  Estimators  of  the 
Spectral  Density.  Soviet  Math.  Dokl.,  Vol.  25,  No.  2,  415-419. 

[2]  Brillinger,  D.  R.  (1981).  Time  Series:  Data  Analysis  and  Theory,  Holden-Day,  San 
Francisco. 

[3]  Burg,  J.  P.  (1967).  Maximum  Entropy  Spectral  Analysis.  Paper  Presented  at  the  37th 
Annual  international  S.E.G.  Meeting,  Oklahoma  City,  Oklahoma. 

[4]  Burg,  J.  P.  (1972).  The  Relationship  between  Maximum  Entropy  Spectra  and  Maxi¬ 
mum  Likelihood  Spectra.  Geophysics,  37,  375-376. 

[5]  Daubechies,  I.  (1992).  Ten  Lectures  on  Wavelets.  SIAM.  Philadelphia. 

[6]  Davis,  H.  T.  (1935).  Tables  of  The  Mathematical  Functions,  Vol  II,  Trinity  University. 

[7]  Donoho,  D.  (1992).  De-Noising  By  Soft-Thresholding.  Technical  Report,  Department 
of  Statistics,  Stanford  University. 

[8]  Donoho,  D.  and  Johnstone,  I  (1992a).  Minimax  Estimation  via  Wavelet  Shrinkage. 
Technical  Report,  Department  of  Statistics,  Stanford  University. 

[9]  Donoho,  D.  and  Johnstone,  I  (1992b).  Adapting  to  Unknown  Smoothness  via  Wavelet 
''  Shrinkage.  Technical  Report,  Department  of  Statistics,  Stanford  University. 

[10]  Donoho,  D.  and  Johnstone,  I  (1992c).  Ideal  Spatial  Adaptation  by  Wavelet  Shrinkage. 
Technical  Report,  Department  of  Statistics,  Stanford  University. 

[11]  Donoho,  D.,  Johnstone,  I.,  Kerkyacharian,  G.  and  Picard,  D.  (1992).  Wavelet  Shrink¬ 
age:  Asymptopia?  Technical  Report,  Department  of  Statistics,  Stanford  University. 


18 


[12]  Gao,  Hong- Ye  (1992).  Spectral  Density  Estimation  by  Using  Wavelets.  Teclinical  Re¬ 
port,  Department  of  Statistics,  Stanford  University. 

[13]  Harvey,  A.  C.  (1981).  Time  Series  Models^  Philip  Allan. 

[14]  Harvey,  A.  C.  (1989).  Forecasting,  Structural  Time  Series  Models  and  the  Kalman 
Filter.  Cambridge  University  Press. 

[15]  Johnstone,  I.,  Kerkyacharian,  G.  and  Picard,  D.  (1992).  Estimation  d’une  densite  de 
probabilite  par  methode  d’ondelettes.  Comptes  Rendus  Acad.  Sciences  Paris  (A),  315 
211-216. 

[16]  Moulin,  P.  (1993a)  A  Wavelet  Regularization  Method  for  Diffuse  Radar- Target  Imaging 
and  Speckle-Noise  Reduction.  Journal  of  Math.  Imaging  and  Vision,  Special  Issue  on 
Wavelet,  Vol.  3,  123-134. 

[17]  Moulin,  P.  (1993b)  Wavelet  Thresholding  Techniques  for  Power  Spectrum  Estimation. 
IEEE  Trans,  on  Signal  Processing,  Vol.  42,  No.  11,  pp.  3126-3136. 

[18]  Parzen,  E.  (1974).  Some  Recent  Advances  in  Time  Series  Modelling.  IEEE  Trans. 
Automatic  Control,  AC-19,  723-729. 

[19]  Priestley,  M.  B.  (1981).  Spectral  Analysis  and  Time  Series,  Academic  Press,  London. 

[20]  Taniguchi,  M.  (1979).  On  Estimation  of  Parameters  of  Gaussian  Stationary  Processes. 
J.  Appl.  Probab.,  16,  575-591. 

[21]  Taniguchi,  M.  (1980).  On  Estimation  of  the  Integrals  of  Certain  Functions  of  Spectral 
Density.  J.  Appl.  Probab.,  17,  73-83. 

[22]  Tukey,  J.  W.  (1978).  When  Should  Which  Spectrum  Approach  Be  Used?  Paper  Pre¬ 
sented  at  the  Institute  of  Statisticians  Conference,  King’s  College,  Cambridge,  July, 
1978. 

[23]  Wahba,  G.  (1980).  Automatic  Smoothing  of  the  Log  Periodogram.  Journal  of  the 
American  Statistical  Association,  75,  122-132. 

[24]  Wittwer,  G.  (1986)  On  the  Distribution  of  the  Periodogram  for  Stationary  Random 
Sequences  (with  Discussion).  Math.  Opera,  und  Stat.,  17,  201-219. 

I 


19 


